#include<iostream>
using namespace std;
#include<cmath>
#include<fstream>
#include"library.h"

double E;

void derivative1(int n, double x, double *y, double *dydx){
	dydx[0] = y[1];
	dydx[1] = - 2 * E * y[0];
}

void derivative2(int n, double x, double *y, double *dydx){
	dydx[0] = y[1];
	dydx[1] = - 2* (E-0.1) * y[0];
}

double psib(double x){

	E=x;

	int N=10000,n=2;
	double a=-1, ya=0, deriv_a=1, b=1;
	int i,j;
	double h = (b-a)/(2*N);
	double *y = new double [n];
	double *ynext = new double [n];
	y[0]=ya; y[1] = deriv_a;
	for(i=0;i<N;i++){
		rk4(a+i*h, h, n, y, ynext, derivative1);
		for(j=0;j<n;j++) y[j] = ynext[j];
	}
	for(i=0;i<N;i++){
		rk4(a+(N+i)*h, h, n, y, ynext, derivative2);
		for(j=0;j<n;j++) y[j] = ynext[j];
	}
	return y[0];
}

int main(){

	double eig = bisection(psib, 0, 2, 1E-15);
	cout<<"Numerical: Eg.s.="<<eig<<endl;
	cout<<"First-order peturbation theory: Eg.s. = "<<M_PI*M_PI/8+0.05<<endl;

	return 0;
}
